Computational reactive–diffusive modeling for stratification and prognosis determination of patients with breast cancer receiving Olaparib

Mathematical models based on partial differential equations (PDEs) can be exploited to handle clinical data with space/time dimensions, e.g. tumor growth challenged by neoadjuvant therapy. A model based on simplified assessment of tumor malignancy and pharmacodynamics efficiency was exercised to discover new metrics of patient prognosis in the OLTRE trial. We tested in a 17-patients cohort affected by early-stage triple negative breast cancer (TNBC) treated with 3 weeks of olaparib, the capability of a PDEs-based reactive–diffusive model of tumor growth to efficiently predict the response to olaparib in terms of SUVmax detected at 18FDG-PET/CT scan, by using specific terms to characterize tumor diffusion and proliferation. Computations were performed with COMSOL Multiphysics. Driving parameters governing the mathematical model were selected with Pearson's correlations. Discrepancies between actual and computed SUVmax values were assessed with Student’s t test and Wilcoxon rank sum test. The correlation between post-olaparib true and computed SUVmax was assessed with Pearson’s r and Spearman’s rho. After defining the proper mathematical assumptions, the nominal drug efficiency (εPD) and tumor malignancy (rc) were computationally evaluated. The former parameter reflected the activity of olaparib on the tumor, while the latter represented the growth rate of metabolic activity as detected by SUVmax. εPD was found to be directly dependent on basal tumor-infiltrating lymphocytes (TILs) and Ki67% and was detectable through proper linear regression functions according to TILs values, while rc was represented by the baseline Ki67-to-TILs ratio. Predicted post-olaparib SUV*max did not significantly differ from original post-olaparib SUVmax in the overall, gBRCA-mutant and gBRCA-wild-type subpopulations (p > 0.05 in all cases), showing strong positive correlation (r = 0.9 and rho = 0.9, p < 0.0001 both). A model of simplified tumor dynamics was exercised to effectively produce an upfront prediction of efficacy of 3-week neoadjuvant olaparib in terms of SUVmax. Prospective evaluation in independent cohorts and correlation of these outcomes with more recognized efficacy endpoints is now warranted for model confirmation and tailoring of escalated/de-escalated therapeutic strategies for early-TNBC patients.

Baseline timepoint TILs Tumor-infiltrating lymphocytes TNBC Triple negative breast cancer Breast cancer (BC) is the most common cancer in women and the most frequent cause of death by cancer in this sex 1 . Fortunately, in the last decades, new and escalated treatment strategies in early-stage disease have led to a substantial reduction in recurrence rates and improvements in survival 2,3 . However, treatment costs and toxicities have also increased substantially 4,5 . This has led to focus new research efforts in better personalizing therapeutic approaches, so to spare unnecessary toxicities and optimize BC care costs, as also advocated by the European society for medical oncology (ESMO) and the broader scientific community in recent years [6][7][8] . In this perspective, the development of tools capable of predicting tumor progression and response to novel therapies might help implementing therapeutic personalization and better identifying patients that might be spared longterm chemotherapy.
In the last few years, mathematical modeling has been entering the arena of oncological research in an attempt to predict spatial and temporal evolution of tumors transferring in-silico models to clinical research and practice 9,10 . Gompertzian and logistic mathematical models were first used to represent tumor cells' growth and invasiveness and have been successively adopted in more sophisticated and complex models for tumor proliferation studies 11,12 . Accumulating evidence is showing that mathematical models based on partial differential equations (PDEs) are potentially exploitable to handle clinical data with spatial dimensions not solely depending on time; which is, for example, the case of tumor growth challenged by neoadjuvant therapy (NAT) 9,[13][14][15] . The PDE approach based on reaction-diffusion models is often employed for cancer modeling. These models define the diffusion and proliferation of the various tumor components, including cancer cells, healthy cells, extracellular matrix etc. with specific mathematical formulas 16 . In this perspective, we preliminarily showed in a restricted cohort of 3 BC patients undergoing NAT, that a computational mass transfer modeling based on a set of PDEs applied at the tumor dynamics might represent a powerful in silico tool to virtualize tumor progression and predict tumor dynamics in response to therapy at the single-patient level 17 .
We hereby retrospectively applied our reactive-diffusive PDEs-based model to a wider BC patient cohort prospectively enrolled in a window-of-opportunity trial at our Institution 18 , to further test whether the combination of personalized diagnostic imaging and clinicopathological tumor/patient variables in mathematical modeling can accurately predict early-stage BC progression and its competition with a suitable NAT for a better patient-adapted planification of the therapeutic strategy.

Methods
Study population. This analysis was retrospectively performed on the BC patients enrolled within the OLTRE "window of opportunity" trial (NCT02681562) with available data for the mathematical modeling. Within the OLTRE study, conducted at the ASST Cremona between 2016 and 2019, treatment-naïve patients with locally advanced non-metastatic HER2-negative BC, with or without a germline BRCA1/2 (gBRCA) mutation, received the PARP inhibitor (PARPi) olaparib at a dose of 300 mg orally for 21 consecutive days, before starting the standard neoadjuvant chemotherapy (CT) 18 . The main objective of the trial was to explore the biological effects of a short course of olaparib, especially in locally advanced triple negative BC (TNBC) independently of the gBRCA status. www.nature.com/scientificreports/ All patients underwent a 18 FDG-PET/CT scan at baseline and after 3 weeks of olaparib ± 3 days; and clinical assessments were conducted at baseline and every 3 weeks ± 3 days. Clinical responses were evaluated through physical exam with caliper and assessed according to RECIST1.1 criteria 18,19 . The same operator performed all physical examinations pre/post olaparib to identify clinical responders (complete response + partial response [CR/PR]) and non-responders (stable disease + progressive disease [SD/PD]). 18 FDG-PET/CT was used to detect radiometabolic responsiveness to olaparib by taking into account baseline and post-olaparib maximum standard uptake value (SUV max ) values for the primary lesion. The same radiologist evaluated all PET/CT responses. Full study details are reported elsewhere 18 . Only TNBC patients were included in the present analysis.
Study hypotheses and objectives. The main hypothesis behind this sub-analysis of the OLTRE trial was to test an in-silico reactive-diffuse model based on PDEs to predict the BC metabolic response to NAT with olaparib alone, as detected by 18 FDG-PET/CT in terms of SUV max after 3 weeks of neoadjuvant olaparib. To differentiate between actual and predicted SUV max values in silico, we will refer to the latter as SUV* max . As a consequence, for the purpose of the present study, tumor dimensional assessments were not considered for the development of our mathematical model. Study procedures: theoretical premises. From the engineering point of view, BC is a single-phase (solid) biomaterial, featuring sharp boundaries delineating the cancer cells population Ø c , that grows and invades a region of interest (ROI) 10 . During NAT, a given mass rate of olaparib Ø d was administered. When Ø c represents an index of metabolic activity, its integration in the ROI can be compared with the measurement of SUV max by 18 FDG-PET/CT and a Gompertzian logistic function can be employed to describe its change rate, as reported in equation Eq. (1), as we and others preliminarily demonstrated 11,20-22 . In the equation, r c represents the nominal personalized biological conversion rate depending on the nanoscale (genomics) and the microscale (cell signals/molecular biology), such as BC invasiveness, aggressiveness or malignancy. In our analysis, r c reflects the growth rate of metabolic activity as detected by SUV max , surrogate of BC malignancy. The 1/r c is a timescale constant, representing the growth rate of Ø c . The variable t stands for time, while the parameter K is an arbitrary constant representing a surrogate of the carrying capacity of the biological matrix. Namely, the limiting nutrients for the onset of the metabolic conversion in the confined space where the BC lesion is detected. The value of K is taken such that the sigmoid function described by Eq. (1) approaches its future asymptote very gradually. Following the Gompertzian function reported in Eq. (1), the tumor growth dynamics can be subdivided into three phases, represented in Fig. 1A. A Phase I of free tumor proliferation (FP), www.nature.com/scientificreports/ when the lesion starts growing from an unknown time (t i ) in the past, until it is first diagnosed at time 0 (t 0 ). The following Phase II is also called challenged proliferation (CP), representing the time when a drug is administered and, if effective, impairs tumor growth/metabolic activity. Usually, a phase III of further tumor growth/increase in metabolic activity is observed anytime a treatment is ceased or resistance is developed. Surgical resection or a treatment change interrupts this phase. Since we aimed at predicting the tumor metabolic activity as detected by 18 FDG-PET/CT in terms of SUV max after 3-week neoadjuvant olaparib alone, we considered as t0 the time of baseline SUV max detection.

Governing equations and model assumptions.
Two PDEs for reaction-diffusion were used in our mathematic model, represented by Eqs. (2) and (3) 20 .
Equation (2) represents the tumoral biomass transport, in terms of dimensionless metabolic activity indicator, while Eq. (3) the drug transport, in terms of olaparib concentration. The source terms R c and R d are defined by the following equations: Importantly, the following assumptions were taken: (1) the physical and functional lag consisting in multiple compartments mediating drug delivery to the tumor lesion was not taken into account, since the effect of mediating compartments would be approximately the same for each single patient and the study had an overall a short time span (3 weeks of olaparib administration) 23 ; (2) a detailed action of the extra-cellular matrix (ECM) was not considered in this model. However, the crosstalk between cancer-associated fibroblasts (CAFs) and tumor cells was taken into account in the bulk effect of the effective diffusion coefficient of tumoral biomass, D c in Eq. (2) 24,25 ; (3) D d in Eq. (3) represented the diffusion coefficients for the drug, i.e. olaparib. For simplicity, we considered the same D values for all patients enrolled into the study (i.e. D c = 1 e−13 m 2 /s; D d = 1e −5 m 2 /s) 26 .
Moreover, to fully understand the above-mentioned equations, the following definitions were adopted: • ε PD : nominal drug efficiency, or aggregated personalized pharmacodynamic (PD) behavior of the drug (i.e. olaparib). It represents the effect of the drug on the tumor in terms of either tumor shrinkage or, in this study, SUV max reduction; • f(t): the indicator of therapy regimen for the administered drug, equal to 1 during the entire neoadjuvant treatment duration. Therefore, the drug concentration in the patient blood was assumed to be constant; • ε PK : the known effect of the clearance, or pharmacokinetic (PK) behavior for the administered olaparib, based on the available drug specifications (i.e. olaparib nominal plasma clearance: 7 l/h 27,28 ). In other words, ε PK brings purposely the PK effect into the PDE representing the distribution of drug concentration. To this end, the drug's nominal plasma clearance (converted in m 3 /s) must be multiplied for the patient's nominal density (in kg/m 3 , conventionally assumed equal to water) then divided for the patient's mass (in kg), in order to obtain ε PK 's desired unit of 1/s, to reach unit consistency in the PDE. This parameter is essential to calculate R d , as previously reported.
Due to the integration of the PDEs system described by Eqs. (2) and (3), the Ø c and Ø d evolve in space and time. The space integration is performed in the available breast volume, while the variables progress over time (t), starting from the in-silico starting time of BC lesion (t i ) to the end of olaparib NAT.
Study procedures: analysis. The system of equations applied to the available breast volume, Eqs. (2-3), supplemented by source term R c and R d definitions equations, along with their initial and boundary conditions previously described, were integrated with the Finite Element Method, by using the COMSOL Multiphysics platform 29 . An unstructured meshing technique was used, yielding for a homogeneous tetrahedral element grid. After a grid independency test, a final mesh of approximately 10.000 elements was employed, to optimize result accuracy and computational times 17 . Direct solver PARDISO was employed as the algebraic engine, while the BDF method was applied for the temporal dependence 17 . Execution durations, for each patient, did not exceed 30 min on a Pentium® Xeon server (Windows® 10 [Microsoft, Redmont, WA, USA], Eightcore-32N at 2.4 GHz, 128 GB RAM) running in serial mode.
Population characteristics were assessed through standard descriptive statistics and variations in mean for main pathological variables pre/post olaparib were assessed through Students' t-test for paired samples. The correlation of the main clinicopathological features (i.e. Ki67%, tumor-infiltrating lymphocytes [TILs] %, age in years, tumor size in mm, PD-L1 TPS and IC %) with the difference between post-and pre-olaparib SUV max values, were used to identify the parameters better correlating with tumor malignancy r c to provide the necessary estimations within the mathematical model. Finally, the post-olaparib SUV max and SUV* max were compared with both parametric (unpaired Student's t-test) and non-parametric (Wilcoxon rank sum test) tests to assess the

Results
The population of the OLTRE trial has been already described elsewhere 18  www.nature.com/scientificreports/ of 17 patients to predict SUV max response. An inspection on the analytical structure of the model revealed that, once applied the known f(t) and ε PK for the specific patient, the progress in the FP phase (Fig. 1A) depended on the previously mentioned t i and r c parameters. Therefore, we firstly performed multiple Pearson correlation analyses to identify the clinicopathological factors potentially associated with SUV max modifications. We noted the SUV max reduction obtained with olaparib was significantly inversely correlated with baseline Ki67 and TILs levels (r = − 0.75 and r = − 0.74 , respectively, p < 0.001 both) only (Fig. 1B). Hence, tumor malignancy, identified with the parameter r c (i.e. SUV max growth rate) was subsequently defined as Ki67/TILs ratio. This mathematical definition was based on the evidences showing in early-stage TNBC that higher values were associated with better survival outcomes [30][31][32] . We normalized Ki67 and TILs values in order to obtain a value range comprised between 0 and 1. Then, the lesion starting time ti was tweaked in each case to closely match each measured baseline SUV max , at the end of the FP phase. Next, for each single patient the models were iteratively run by manually-adjusting the nominal personalized olaparib efficiency (ε PD ) in order to come up, in the subsequent CP phase, with minimal relative deviations in the computed post-olaparib SUV* max volumes, with respect to the corresponding actual post-olaparib SUV max measurements. Results are summarized in Table 2.
At this point, we had a very good reproduction of each clinical case, but no actual potential of prospective evaluation, yet. We then plotted ε PD against r c (i.e. Ki67/TILs ratio) discovering that two sub-cohorts appeared neatly depending on the range of normalized TILs level (Fig. 1C). The patient normalized Ki67 and TILs values could be grouped along two different interpolating curves based on the following functions: (1) for TILs level below the value of 0.5: ε PD = 0.1227 · r c 1.5184 , with R 2 = 0.9277 (2) for TILs level above the value of 0.5: ε PD = 0.004 · r c 2 + 0.0105 · r c + 0.0252, with R 2 = 0.9448 Hence, with ε PD expressed with the above dependences, all of the driving parameters in Eqs. (2-3) were identified, and the model could be solved in a full prospective mode, with its solution depending on the available baseline SUV max and baseline Ki67 and TILs levels. Consequently, in a second computational stage the models were iteratively run again, by keeping the same optimized values for t i and r c reported in Table 1. The computed SUV* max , along with the relative deviations with respect to the corresponding true SUV max measured by PET, are listed in Table 3.
The prediction quality was adequate, with computed post-olaparib SUV* max differing of < 1 unit with respect to the actual post-olaparib SUV max for 12 (70.6%) patients and in the range of 1-2 units for the remaining 5 (29.4%), independently from gBRCA status. The correlation between post-olaparib SUV max and computed SUV* max was positive and very high, according to parametric and non-parametric methods, as well (Pearson's r = 0.95, p < 0.0001, Spearman's rho = 0.93, p < 0.0001). Importantly, the numerical difference between postolaparib SUV max and SUV* max was not statistically significant by using both parametric (p = 0.813, p = 0.866 and p = 0.856) and non-parametric (p = 0.945, p = 0.841 and p = 0.977) statistics for the overall, gBRCA-mutant and gBRCA-wild-type populations, respectively.
Finally, we exercised the model in a virtual scenario of different pharmacodynamic efficiency (ε PD ) and olaparib duration, for a randomly selected patient of our cohort (patient n.13), with the objective of bringing Table 2. Detailed values of the driving parameters of the mathematical model, with clinical and computed SUV max . t i , lesion starting time; r c , biological conversion rate. In this case each r c is the results of Ki67/TILs; ε PD , olaparib efficiency (effect of the drug on the body); SUV max , actual maximum standard uptake volume measured by 18  www.nature.com/scientificreports/ out the non-linear relationship between the predicted outcome (SUV* max ) and time. As observable in Fig. 2, the progress of the tumor computed SUV* max for patient n.13 was compared, with the variable values leading to the results reported in Table 3 (Case A), to a virtual case represented by the same patient with a 20% decrement of ε PD and an increase in 6 days of olaparib duration (Case B). We observed that in case of variation of both response to NAT and NAT duration, it was again possible to achieve the same quantitative reduction in predicted SUV max (i.e. SUV* max = 1.72) by stretching the NAT duration (Fig. 2). Such results were obtained by executing the model iteratively, until the same original volume of Case A was achieved. It is confirmed therefore the results obtainable from the proposed model, through its non-trivial solution, are strongly intertwined and non-linear in nature, and the model presents the flexibility and adaptability to potentially tailor in silico an entire therapeutic strategy at the patient level. Table 3. Values of εPD for each patient, with related SUV* max and their absolute deviations with respect to the clinical measurements of post-olaparib SUV max . SUV* max , predicted SUV max values; εPD, olaparib efficiency (effect of the drug on the body); SUV, standard uptake volume.

Discussion
We tested, in a cohort of 17 patients affected by early-stage TNBC treated with 3 weeks of olaparib in a "window of opportunity" trial, a mass transfer PDEs-based reactive-diffusive mathematical model of tumor growth which might be capable of efficiently predicting the response to olaparib in terms of SUV max detected at 18 FDG-PET/ CT scan. The model showed, without any preliminary assumption, the effective pharmacodynamic efficiency of olaparib was strongly dependent on basal TILs level and SUV max growth rate. The latter was represented by a mathematical parameter that in our case was directly dependent on Ki67 expression and TILs count. By knowing the basal SUV max , Ki67% and TILs levels it would be possible with this approach to predict with a very small margin of error the SUV max change after 3 weeks of olaparib. Although, this is not standard of care, it opens up to the possibility of early-detecting a drug efficacy (olaparib in this case) and investigating the personalization of escalated or de-escalated therapeutic strategies in early-stage TNBC in future research.
The need for biomarkers to guide escalated/de-escalated and more personalized therapeutic approaches is a hot topic in current Oncology 33,34 . Moreover, TNBC are historically a subgroup of breast tumors with poor prognosis and lack of biomarkers for effective personalized therapeutic approaches 35,36 . The most notable and very recent (for their therapeutic implications) exceptions are the determination of gBRCA status, in both adjuvant and metastatic setting, to decide whether or not to prescribe PARPi and the assessment of PD-L1 levels/positivity for the prescription of 1 st -line immunotherapy + chemotherapy in metastatic disease [37][38][39] . However, the main issues with biomarker discovery research include: (1) higher difficulty related to lack of funding or difficult access to data from clinical trials, (2) regulators are less prone/used to approve prognostic and predictive biomarkers than novel therapeutic options, (3) the high number of false discoveries, (4) lack of reproducibility or complex/ inviable implementation in clinical practice [40][41][42] . This translates into a slow and defective transfer of knowledge from the laboratory to the clinical research and/or practice scenario.
Noteworthy, the evaluation of Ki67 is an already established prognostic biomarker in BC, which expression is strongly associated to tumor proliferation and growth 43 . Higher levels are usually associated to higher tumor aggressiveness and worse prognosis, although in TNBC there is no established cut-off to define high vs. low Ki67 levels, differently from hormone receptor-positive BC 35,43,44 . Nevertheless, a recent meta-analysis of 35 independent studies (~ 8000 patients with resected TNBC) suggested that a cut-off of 40% would be associated with higher recurrence risk and mortality 45 .
The morphological evaluation of TILs in BC has gained attention in the last few years, when preliminary evidences started to show a potential prognostic and predictive role, especially in TNBC and HER2-positive BC 46 . A recent retrospective analysis on more than 2000 patients showed a clear favorable prognostic role for higher TILs levels (≥ 30%) in early-stage TNBC, independently from main clinicopathological factors 30 . This evidence adds to the strong independent prognostic role showed by TILs in residual disease after neoadjuvant chemotherapy 47 and a retrospective analysis where higher TILs were found to be independently associated to multiple survival endpoints in patients from an old cohort not treated with (neo)adjuvant chemotherapy. In the same study, stage I tumors with TILs ≥ 30% showed a 5-year overall survival of 98% 48 . These evidences assure our mathematical model was built on biologically and clinically meaningful parameters. Furthermore, our model might be a powerful tool for the personalization of BC care not necessarily requiring the detection of novel costly biomarkers.
There are several limitations that will have to be overcome to promote the implementation of this mathematical framework in the clinical research and practice scenarios. First of all, we had the possibility to study our mathematical model in a cohort of 17 patients, which is more than what is usually done in this research field, where modelling frameworks based on single patients are the norm 13,17,21 . However, wider cohorts are required to validate our findings. Secondly, olaparib is still not approved in the neoadjuvant setting, nor in gBRCA-wild type TNBC. Nevertheless, PARPi showed activity and efficacy in several solid tumors also beyond gBRCA mutational status 38 and the model performed quite well in both gBRCA-wild-type and mutant patients. Moreover, another PARPi, namely talazoparib, already showed promising neoadjuvant efficacy in gBRCA-mutant TNBC 49 and the same olaparib is currently under evaluation in a phase II study in monotherapy or in combination with the immune-checkpoint inhibitor durvalumab in early-stage HER2-negative BC with either a germline or somatic BRCA1/2 mutations (OlimpiaN, ClinicalTrials.gov identifier: NCT05498155). In addition, also the PARPi niraparib recently showed in a pilot study, a high tumor response (90.5%) on MRI along with high pCR rate (40.0%) in gBRCA-mutant HER2-negative breast cancer patients who received it as monotherapy in the neoadjuvant setting 50 . This means that PARPi have the potential to become a neoadjuvant therapeutic option for TNBC in the next future and the mathematical model hereby tested might be envisioned as an in silico predictor of response for further upfront implementation of escalated (e.g. the addition of immunotherapy and/or addition of posterior chemotherapy) or de-escalated strategies (e.g. PARPi monotherapy alone). Yet, this should be extensively tested in the future. For the present, it will be interesting to understand if the same model can be applied to different available therapies, for example to evaluate the opportunity to add carboplatin to standard anthracycline-taxane-based neoadjuvant chemotherapy, to avoid the cardiotoxic anthracyclines or even the evaluate the addition of immune-checkpoint inhibitor pembrolizumab. A better correlation with clinical outcomes has to be further tested, as well.
Another limitation, is that SUV max detected by 18 FDG-PET/CT is not the standard of care for assessing response to neoadjuvant therapy in BC. However, recent results of the PHERGain trial in early-stage HER2 + BC, showed that ~ 1/3 of patients with HER2 + BC treated with neoadjuvant trastuzumab and pertuzumab might be spared chemotherapy, thanks to the degree of metabolic response detected with FDG-PET/CT after just 2 cycles of the anti-HER2 combination and subsequent type of pathologic response detected after surgery 51 . In this line, our study has to be intended as proof-of-concept analysis where our aim was to essentially prove that a PDEs reactive-diffusive mathematical model based on the assumption of Gompertzian growth could be applied not www.nature.com/scientificreports/ only to a setting where tumor dimension/growth is directly considered, but also to demonstrate the capability to track tumor metabolism changes as detected by a standardized and objective measurement parameter (i.e. 18 FDG-PET SUV max ) as a surrogate of malignancy and relate it to tumor response to treatment. In this perspective, we believe we succeeded in our intent, though replicating our findings in other independent but similar databases will be crucial. Finally, the software COMSOL Multiphysics is not open-access and a more user-friendly interface should be envisioned for a broader implementation outside a pure Engineer/Mathematic environment. Nevertheless, COMSOL's computational robustness is already consolidated and the proposed mathematical model has the advantage of potentially being run at low cost on standard desktop computers, being also virtually adaptable to any proliferation/therapy scenarios, as also preliminarily observed in Lymphomas and a different BC setting 17,21 .
To conclude, we observed that a mathematical framework based on realistic multidimensional governing PDEs, could be directly applied to the tailored simulation of an early therapeutic response to the PARPi olaparib in early-stage TNBC by using 18 FDG-PET/CT scan, at the single patient level. The analytical and computational structure of the model sets the basis for further development in in-silico prognosis in Oncology, where the model has the potential to be tested in any virtual scenario on any possible patient, for any combination of the variable space in a sustainable way, to inform and support the Oncologists in their therapeutic decisions.
Prospective evaluation in independent cohorts and correlation of these outcomes with more recognized efficacy endpoints is now warranted.

Data availability
The datasets generated during and/or analyzed during the current study are available from the Corresponding Author upon reasonable request.